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1 Linear stability analysis 

Equilibria are not always stable. Since stable and unstable equilibria play quite 
different roles in the dynamics of a system, it is useful to be able to classify equi¬ 
librium points based on their stability. 

Suppose that we have a set of autonomous ordinary differential equations, 
written in vector form: 

x = f(x). (1) 

Suppose that x* is an equilibrium point. By definition, f(x*) = 0. Now sup¬ 
pose that we take a multivariate Taylor expansion of the right-hand side of our 
differential equation: 
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The partial derivative in the above equation is to be interpreted as the Jacobian 
matrix. If the components of the state vector x are (xi,X2, • • • ,x n ) and the compo¬ 
nents of the rate vector f are (/i, /2, • • •, /«), then the Jacobian is 
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Now define 8x = x — x*. Taking a derivative of this definition, we get 8x = x. If 
8x is small, then only the first term in equation 2 is significant since the higher 
terms involve powers of our small displacement from equilibrium. If we want to 
know how trajectories behave near the equilibrium point, e.g. whether they move 
toward or away from the equilibrium point, it should therefore be good enough to 
keep just this term. 1 Then we have 

8x = J* 8x, 

where J* is the Jacobian evaluated at the equilibrium point. The matrix J* is a 
constant, so this is just a linear differential equation. According to the theory 
of linear differential equations, the solution can be written as a superposition of 
terms of the form where { Xj } is the set of eigenvalues of the Jacobian. 

The eigenvalues of the Jacobian are, in general, complex numbers. Let Xj = 
Hj + iVj, where pj and Vj are, respectively, the real and imaginary parts of the 
eigenvalue. Each of the exponential terms in the expansion can therefore be writ¬ 
ten 

The complex exponential in turn can be written 

e ,Vjt — cos (Vjt) +ism(Vjt). 

The complex part of the eigenvalue therefore only contributes an oscillatory com¬ 
ponent to the solution. It’s the real part that matters: If pj > 0 for any j, e^d 
grows with time, which means that trajectories will tend to move away from the 
equilibrium point. This leads us to a very important theorem: 

Theorem 1 An equilibrium point x* of the differential equation 1 is stable if all 
the eigenvalues of J*, the Jacobian evaluated at x*, have negative real parts. The 
equilibrium point is unstable if at least one of the eigenvalues has a positive real 
part. 

Because we are only keeping a locally linear approximation to the vector field, an 
analysis based on this theorem is called a linear stability analysis. 

Note that the theorem is silent on the issue of what happens if some of the 
eigenvalues have zero real parts while the others are all negative. This case can’t 
be decided based on linear stability analysis. The nonlinear terms we left out 

'This assumes that the Jacobian evaluated at the equilibrium point isn’t, in a sense to be deli ned 
later, zero. 
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of equation 2 in fact determine the stability in this case. Dealing with this case 
requires a nonlinear theory which we discuss later. 


Example 1.1 Let us return to the Lindemann mechanism, for which phase-plane 
analysis has already shown us that the equilibrium point is stable. The di¬ 
mensionless ODEs are 


a — —a~ + aab, 
b — a — a ab — b. 


The equilibrium point is (0,0). The Jacobian matrix is 


d a 
da 

d (1 
db 


—2a + ab 

a a 

db 

- da 

db 
db - 


2 a — ab 

—a a — 1 


Evaluating the Jacobian at the equilibrium point, we get 



The eigenvalues of a 2 x 2 matrix are easy to calculate by hand: They are 
the solutions of the determinant equation 


|AI-J| =0. 


In this case, 


X 0 
0 X+l 


= X(k+l) =0. 


The solutions of this equation can be read by inspection: A, = 0orA, = —1. 
One of the eigenvalues is zero, so we can’t tell from the linear stability 
analysis alone whether or not the equilibrium point is stable. Of course, we 
already know from the phase-plane analysis that it is. 


Example 1.2 The law of microscopic reversibility says that we can’t have truly 
irreversible elementary chemical reactions, although this might be a good 
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approximation if the reaction is strongly product-favored. Consider there¬ 
fore the fully reversible Lindemann mechanism: 

k\ 

A + A A + B, 

k~ 1 

h 

B ^ P. 

k-2 

Recall that A + B + P^Aoisa constant in this mechanism. Defining 

a — k\A/k2 , CL = k-i/k\, 

b = k\Bjk 2l $ = k- 2 /k2, 

x = &2C ao=kiAo/k2, 

and using the conservation relation, the dimensionless ODEs become 

a — —a~ + aab, 
b — a 2 — aab — b + P(«o — a — b). 


If we set a = b = 0, we discover that there are two equilibrium points: 


= ( 0 , 

W) = 


P«o 
1 + P7’ 

afiao 


(3«o 


1 +13( 1 + oc) 1 + p( 1 + oc) 


The Jacobian of the differential equations is 


—2a + ub a a 

2 a — ab — (I —a a — 1 — p 


The Jacobian evaluated at the first equilibrium point is 


ab f 0 
aZ^-p -1 - P 


We could substitute the value of b T directly into J T , but I usually leave this 
for a later step, and then only if it’s necessary to determine the sign of the 
eigenvalues. 
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The eigenvalues satisfy the equation 

Hccfct - P ) ^+1 + P =( X - a ^ t )( ?l + 1 + P ) = a 

The eigenvalues associated with (a',£>') are therefore 

Xi = a b\ 
and ^2 = —(1 + |3). 

Since b ^ is positive, Xi is clearly positive. It follows that is an 

unstable equilibrium point. 

We could try to work out the stability of the other point by hand, but it’s 
messy. In this case, it’s far better to use Maple. The steps in the analysis 
are much the same, although it takes a few tricks to get to the bottom of this 
exercise. We start by defining the differential equations: 

> adot := (a,b) -> -a~2+alpha*a*b; 

adot := (a, b) —> — a 2 + aab 

> bdot := (a,b) -> a~2 - alpha*a*b - b + beta*(aO-a-b); 

bdot := ( a , b) —> a 2 — aab — b + P (aO — a — b ) 

Define the steady state: 



> simplify(adot(astar,bstar)) ; 
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0 


> simplify(bdot(astar,bstar)) ; 

0 

Calculate the Jacobian: 


> with(linalg) : 

Warning, the protected names norm and trace have been 
redefined and unprotected 

> J := jacobian ([adot (a,b), bdot (a,b) ], [a, b]) ; 


J:= 


—2a + <Xb a a 

2a — ab — p — aa— 1 — (3 


Substitute the steady state into the Jacobian: 


> Jstar := map(x->subs(a=astar,b=bstar, x) , J) ; 


afiaO arfiaO 


Jstar := 


%i 


%1 


OtfiaO a 2 fiaO 

- %i %T - 

%1 := l + p(l + a) 

Determine the characteristic polynomial (the polynomial whose roots are 
the eigenvalues). Collect the result in powers of the eigenvalue X: 


> collect(charpoly(Jstar,lambda),lambda) ; 

2 (a 2 pr/O + aPr/0+|3a + |3 2 a + |3 2 + 2p + 1)X 
~ + 1+p+pa 

a 2 p 2 aO + a p 2 aO + a p aO 
+ 1+p+pa 

This polynomial is in the form 


At + ciX + cq. 
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Note that the coefficients co and c i are both positive. The roots of this 
polynomial are 

^ = \ |-ci ± \j Cj — 4c 0 |. 

Since co is positive, the quantity under the square root is either smaller 
than C|, or it is negative. If negative, the solutions are complex with real 
part —ci, which is negative. Otherwise, the square root must be smaller in 
absolute value than ci, so that the two eigenvalues must still be negative. 
Either way, we conclude that the steady state is stable since the real parts of 
both eigenvalues must be negative. 

Note that the possibility of a zero eigenvalue disappeared in the last ex¬ 
ample when we considered the effect of the law of microscopic reversibil¬ 
ity. The zero eigenvalue was therefore an artifact of having an incomplete 
chemical model since this feature disappears when we include a nonzero 

value for the rate constant k . 2 , no matter how small this value is. This 

feature of the irreversible Lindemann mechanism is not structurally sta¬ 
ble since an arbitrarily small change to the model makes it disappear. This 
is a minor structural instability since the qualitative behavior of the model 
(approach to the equilibrium point) is pretty much the same, whether the 
leading eigenvalue is zero or negative. However, there are models which 
have severe structural instabilities, i.e. models in which the whole behavior 
of the model is significantly different depending on whether one includes a 
small term or not. Structurally unstable models are generally considered to 
be bad models since we generally don’t know the exact evolution equations 
for a model. If the model’s behavior changes when we make small changes 
to the equations, then we probably can’t trust the model’s predictions very 
well. 

Note also that we have our first example of a model with two equilibrium 
points, one of which is stable while the other is unstable. It is instructive 
to look at the phase portrait, which we’ll draw using Maple this time since 
we have it running anyway. We first have to pick some values for our pa¬ 
rameters. Since our analysis indicates that the stability properties of the two 
fixed points never change, it doesn’t much matter what we pick. 

> alpha := 1; 


a := 1 
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> beta := 1; 


P:=l 

> aO := 1; 


aO 1 

The coordinates of the stable equilibrium point are 

> astar; bstar; 


1 

3 

1 

3 

The unstable equilibrium sits on the b axis at coordinate 

> bdagger := beta*aO/(1+beta); 

bdagger := ^ 

We’re now ready to create the phase portrait. 

> with(DEtools) : 

Warning, the name adjoint has been redefined 

> with(plots): 

Warning, the name changecoords has been redefined 

> pp := phaseportrait([diff(a(t),t)=adot(a(t),b(t)), 
diff(b(t),t)=bdot(a(t),b(t))], [a(t),b(t)], t=0..10, 

[[a(0)=i/5,b(0)=0]$i=0..5, [a(0)=i/5,b(0)=1]$i=0..5, 

[a(0)=0.02,b(0)=0], [a(0)=0.02,b(0)=1]], a(t)=0..a0, 
b(t)=0..a0, arrows=NONE, stepsize=0.01, 
linecolor=green): 

We’re storing the plot in the variable pp. We’ll also plot the nullclines, and 

store them temporarily: 
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> a_nulll := plot(a/alpha,a=0..1,color=red): 

> a_null2 := plot([0,b,b=0..1],color=red,thickness=3): 

> b_null := plot((a~2+beta*(aO-a))/(alpha*a+l+beta), 
a=0..1): 

The display () function puts all four pictures out on one graph: 

> display({pp,a_nulll,a_null2,b_null}, labels=[a,b], 
scaling=CONSTRAINED); 



Note that a — 0 on the b axis. Accordingly, if we start a system with a = 0, 
it moves along the b axis toward the unstable fixed point. However, if we 
start anywhere off the b axis, we eventually end up at the stable fixed point. 
We say that (0, b^) is a saddle point, i.e. a point with one stable and one un¬ 
stable direction. The stable direction corresponds to the negative eigenvalue 
while the unstable direction corresponds to the positive eigenvalue. 

The stable equilibrium point (a*.If ) on the other hand is called a stable 
node. It has two real, negative eigenvalues. 2 

2 We’ve shown that the two eigenvalues always have negative real parts. Proving that they are 
also real is a pain, but it is possible. 
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When the eigenvalues of a fixed point are complex, the point is called a focus. 
Fixed points of planar systems can be almost completely classified based on their 
eigenvalues: 


Eigenvalues 

Fixed point 

Flow 

complex with positive real parts 

unstable focus 


complex with negative real parts 

stable focus 

■s 

\ k 

real and positive 

unstable node 

% 

real and negative 

stable node 


one positive and one negative 

saddle point 

“Av 


The only problem is with systems which have one or more eigenvalues with real 
part zero. In these cases, all of the above behaviors are possibilities, along with 
one other, namely a centre, which is a closed orbit (e.g. a circle) which is neutrally 
stable. The latter means that if we displace the system off its orbit, it will continue 
to circle the equilibrium point, but at a new radius determined by the disturbance. 

Linear stability analysis can be carried out for higher-dimensional systems 
although, predictably, it gets harder to do things analytically. The main ideas 
remain intact however. 


2 Liapunov functions 

Linear stability analysis tells us how a system behaves near an equilibrium point. 
It does not however tell us anything about what happens farther away from equilib¬ 
rium. Phase-plane analysis combined with linear stability analysis can generally 
give us a full picture of the dynamics, but things become much more difficult in 
higher-dimensional spaces. In this section, we consider a technique due to Lia¬ 
punov 3 which can be used to determine the stability of an equilibrium point “in 
the large”, i.e. both near and far from the equilibrium point. 

Liapunov’s method is based on a simple idea. Suppose that V (x) is a function 
of our state variables which has a minimum at an equilibrium point and which has 

? Also frequently spelled “Lyapunov”. 
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no local minima. (Think of a paraboloid.) Now suppose that we can show that 
the dynamics of our system results in a steady decrease in V in some (possibly 
large) neighborhood of the equilibrium point. This necessarily means that we are 
tending toward the minimum of V, which is just the equilibrium point. Having 
shown this, we can conclude that the equilibrium point is stable over the entire 
neighborhood of x* over which V decreases. A function V with these properties 
is called a Liapunov function. 

Let us formali z e this idea: 

Definition 1 Let U be a region of phase space containing the equilibrium point 
x*. Let V : U —> R be a continuous and differentiable function. V is a positive 
definite function for the point x* if it satisfies the following two conditions: 

1. V(x*) = 0 , and 

2. V (x) > 0 for x E 1/ — {x*}. 

Theorem 2 Let x* be an equilibrium point of the differential equation 1, and let 
V be a positive definite function for this point. The equilibrium point is asymptot¬ 
ically stable (the solutions tend to this point 4 ) for initial conditions in the neigh¬ 
borhood U ofx* ifV (x) < 0 for allxeU — {x*}. 

In order to use this theorem, we have to obtain a Liapunov function. Unfortu¬ 
nately, it’s often really difficult to come up with a Liapunov function for a given 
system, except in some special cases where the physics of the problem suggests a 
particular choice. 


Example 2.1 A Hooke’s law spring is subject to a restoring force —kx, where 
x is the displacement from equilibrium. If the spring is pulling around a 
mass m on a lubricated surface, the damping (frictional) force is typically 
proportional to the velocity and opposite in direction. Thus, 

F = ma — —kx — pv. 

4 The technical distinction being made here is between stability, which just means that the solu¬ 
tions don’t wander away from the point, and asymptotic stability, which means that they actually 
tend toward the point with time. 
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By definition, 


v 

and a 


dx 

dt ’ 
dv 

dt 


We therefore obtain the following planar system: 5 


dx 

dt 

dv 

dt 


v, 

— (—kx — in ’). 
m 


The equilibrium point of this system is clearly (x*,v*) = (0,0). In mechan¬ 
ical systems with friction, the total energy is often a Liapunov function. 
Since the force and potential energy 6 V are related by F — —dV/dx, the 
potential energy associated with the Hooke’s law force is jkx 2 . The kinetic 
energy is jmv 2 , so the total energy is 


E 


-kx 2 + -/nv . 
2 2 


The total energy is a positive definite function (zero at the equilibrium point, 
and increasing away fromt this point). Now consider 

g dE dx dE dv 
dx dt dv dt 

— kxv + mv — (—kx — /uv) 
m 

= -/jv 2 < 0 . 

5 Noting that a = d 2 x/dt 2 , this system can also be written as a single second-order ODE 

d 2 x dx 

and this is in fact the approach taken in many mechanics courses. What we are seeing here is 
that second-order ODEs (and higher-order equations, for that matter) are equivalent to sets of 
fi rst-order ODEs, provided we use some of the derivatives (e.g. v = dx/dt) as variables. 

6 The potential energy V should not be confused with the Liapunov function V. The traditional 
notations of mechanics and of Liapunov functions are in conflict here. 
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Since E decreases everywhere except at the equilibrium point itself, the 
equilibrium point is globally asymptotically stable. In order words, starting 
from any initial conditions, we will eventually reach the equilibrium point. 
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